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Abstract 


Investigations of the phase diagram of biaxial liquid crystal systems 
through analyses of general Hamiltonian models within the simplifica¬ 
tions of mean-field theory (MFT), as well as by computer simulations 
based on microscopic models, are directed towards an appreciation 
of the role of the underlying molecular-level interactions to facilitate 
its spontaneous condensation into a nematic phase with biaxial sym¬ 
metry. Continuing experimental challenges in realising such a system 
unambiguously, despite encouraging predictions from MFT for exam¬ 
ple, are requiring more versatile simulational methodologies capable of 
providing insights into possible hindering barriers within the system, 
typically gleaned through its free energy dependences on relevant ob¬ 
servables as the system is driven through the transitions. The recent 
brief report from this group [B. Kamala Latha, et. al., Phys. Rev. 
E 89, 050501(R), 2014], summarising the outcome of detailed Monte 
Carlo simulations carried out employing entropic sampling technique, 
suggested a qualitative modification of the MFT phase diagram as 
the Hamiltonian is asymptotically driven towards the so-called partly- 
repulsive regions. It was argued that the degree of the (cross) coupling 
between the uniaxial and biaxial tensor components of neighbouring 
molecules plays a crucial role in facilitating, or otherwise, a ready con¬ 
densation of the biaxial phase, suggesting that this could be a plau¬ 
sible factor in explaining the experimental difficulties. In this paper, 
we elaborate this point further, providing additional evidences from 
curious variations of free-energy profiles with respect to the relevant 
orientational order parameters, at different temperatures bracketing 
the phase transitions. 
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1 Introduction 

The thermotropic biaxial nematic phase which was predicted by Freiser pQ 
nearly fonr decades ago has attracted considerable attention recently for var- 
ions reasons, ranging from a fnndamental question of conducive experimental 
conditions for its realization to its envisaged applications in display devices. 
Though predictions made by various mean-field (MF) theoretic treatments 
0 - 0. Landau free energy based analyses HDl - H and computer sim¬ 
ulations Ha - 1221 support the feasibility of such a phase, success on the 
experimental front has been rather modest |23]. Experimentally the biaxial 
phase was hrst obtained in a lyotropic, ternary mixture of potassium laurate, 
1-Decanol and D 2 O in 1980 |23] and more recently in bent-core compounds 
PU [2Sj , organo-siloxane tetrapodes EZIEH], LC polymers [29] and colloidal 
systems of Goethite particles I3ni. Though recent experiments [31] 132] point 
to low transition enthalpies for rod-disc systems, an unambiguous biaxial 
phase has not been established in such systems. From the point of view 
of application, it is anticipated that the minor director could switch more 
readily compared to the major director in an external field [231 El] lead¬ 
ing to faster response times. Even in the case of recent bent-core molecules 
there appears to be a debate on the consistency in the experimental Endings 
[351 l36l 1371 l38] . Achieving spontaneous macroscopic biaxiality in nematic 
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liquid crystal phases with appreciable biaxial order appears at the moment 
to be a challenge. 

The recent theoretical studies on the other hand point to a more opti¬ 
mistic picture: they predict that the condensation of a biaxial phase could 
occur over a wide range of the Hamiltonian parameter space of a general 
quadrupolar model 0 - 0. However, the analysis of the mean held model 
was noted to be unsatisfactory, as the phase behaviour of the biaxial system 
in the limit of vanishing intermolecular biaxial interaction traversing in the 
process the so-called partly repulsive region of the Hamiltonian was found to 
be contravening the biaxial phase stability criterion [9]. In this context we 
revisit the mean-held phase diagram with detailed Monte Carlo simulations. 
Main results of this study were briehy presented recently |39]. The other 
MC work on the so-called /x-model [10] was also similarly concerned with 
the consequences of the contribution of a repulsive interaction term in the 
Hamiltonian. 

In this paper, we present the details of a qualitatively diherent type of 
Monte Carlo sampling that we adopted for the study. It was observed that 
the sampling methods to extract equilibrium averages based on equilibrium 
ensembles (constructed using the Metropolis algorithm |TT]) largely lead to 
results in accord with the MFT in the so-called attractive region of the Hamil¬ 
tonian parameter space. Keeping this in view, we adopted the Wang-Landau 
sampling procedure [12] augmented by frontier sampling to determine 

the representative density of states of the system, enabling the calculation 
of all relevant thermodynamic properties. We hnd that this more versatile 
and efficient technique results in qualitatively different results in certain re¬ 
gions of the parameter space, leading to the proposal of a modihed phase 
diagram (relative to MFT). We argue that, such differences, which develop 
progressively as the ’partly repulsive region’ is reached, are important in un¬ 
derstanding the relative roles of different contributions to the intermolecular 
tensor interactions. 

The mean held Hamiltonian model employed and its representation for 
purposes of simulation are outlined in section H. The sampling technique and 
the simulation details are presented in section HI. The observations based 
on these computations are presented in section IV, followed by conclusions 
in section V. 


2 Hamiltonian model 

The MF analysis [3] - [9], is based on the general quadrupolar orientational 
Hamiltonian, proposed by Straley j2] and set in terms of tensors [5]. Ac- 
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cordingly, the interacting biaxial molecules are represented by two pairs of 
symmetric, traceless tensors {q, b) and {q , b ). Here q and q are uniaxial 
components about the unit molecular vectors m and m , whereas b and b 
(orthogonal to q and q , respectively), are biaxial. These irreducible compo¬ 
nents of the anisotropic parts of susceptibility tensor are represented in its 
eigen frame (e,ej_,m) as 


q ■ 

b : 


m® m — 


I 
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e (g) e — e_i_ 0 e_i_ 


(la) 

(lb) 


where I is the identity tensor. Similar representations hold for q' and b' in 
the eigen frame (e , e_|_, m ). The interaction energy is written as 


H = -U[iq- q'+ ^{q-b' + q' ■b) + Xb- b'] (2) 

where U is the scale of energy, ^ = ± 1 , 7 and A are dimensionless interac¬ 
tion parameters, determining the relative importance of the uniaxial-biaxial 
coupling and biaxial-biaxial coupling interactions between the molecules, re¬ 
spectively. 

Mean-field analysis of the Hamiltonian identifies a triangular region OIV 
in the ( 7 , A) plane - called the essential triangle - representing the domain 
of stability into which any physical system represented by Eqn. (|^ can be 
mapped mi (see Fig. [^. The dispersion parabola A = 7^ 0 traverses 
through the interior of the triangle, intersecting IV at the point T, called 
the Landau point. Region of the triangle above the parabola corresponds 
to a Hamiltonian where all the terms are attractive, while the region below 
is noted to be partly repulsive particular, a mean-field (MF) phase 

diagram was predicted i as a function of the arc length OIV (Fig. [^, 
denoted by A*, defined as A* = A on the segment 01, and 

„ (1 + VI37) 

^ - 3 ' 

with 

. (1 - 3A) 

^ 2 

covering the segment IV. The MF phase diagram predicts for A* < 0.22 
(7 = 0, A < 0 . 22 ) a two stage transition from the isotropic to a biaxial phase, 
with an intervening uniaxial nematic phase. The uniaxial-biaxial transition 
is computed to be second order {Nb = Njj — I) upto the point Oi (7 = 0, A ~ 
0 . 2 ), and then changes to hrst order {Nb — Njj — I), till O2 (7 = 0 , A ~ 0 . 22 ). 
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Figure 1: (color online) Essential triangle : Region of biaxial stability. 01 
and IV are uniaxial torque lines intersecting at the point I (0,1/3). OT is the 
dispersion parabola which meets the line IV at the Landau point T. Point V 
(1/2, 0) is the limit of biaxial stability for the interaction. Ci (0,0.2) and 
(5/29,19/87) are tricritical points and C 2 (0.22,0) is a triple point (G. De 
Matteis et al, Continuum Mech. Thermodyn. 19, 1-23 (2007)). K (0.2, 0.2) 
is a point where p = -1 (refer to text). 
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For the rest of the range of A*, a direct isotropic-biaxial transition is expected, 
extending upto V in Fig. This transition is predicted to be hrst order 
{Nb - I) for A* < 0.54 (Cs, A* = 0.54, (7 = 5/29, A = 19/87)) and second 
order {Nb = I) npto the point V (7 = 0.5, A = 0.0). Hence C\ and are 
tricritical points and C 2 is a triple point. 

We consider here the diagonal form of the interaction Hamiltonian in 
Eqn. @ mu expanded as a snperposition of two qnadratic terms, i.e 

H = —U{a'^q~^ ■ + a~q~ ■ q~') 

where q^ and q~ are orthogonal molecnlar biaxial tensors represented as 

= g + 7±h 


with 

, 3A - 1 ± W(3A - 1)2 + 1272 

7 = -- 

67 


and 


7 -7 

7- _^+ 


7-7^ 

7" -7^ 


Along 01 where 7 = 0, = g, = h, a"*" = 1 and a~ = X implying that 

q~^ is pure uniaxial and q~ is pure biaxial and the Hamiltonian reduces to 
an interaction in terms of a single parameter A 


H = -U{q-q' + Xb-b'). 


(3) 


Similarly, along IV, dehned by 1 — 3A — 27 = 0, the Hamiltonian is 
expressed in terms of uniaxial tensor q^ and biaxial tensor as [9] 


H = ■ bi') ( 4 ) 

where 

Q2 = -^Q =(e®e--) 

3 

^2 = -q^ = {fri ® m — ej_ ® ej_) 

and 

(1 - 9A) 

The pair-wise interaction in Eqn. Q now reduces to 
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H = 


e- -) 


" - 3) - 


ej_ — m 0 m) • (e 


e I — m 


m )]. 


(5) 

with U' = f/(l — A)/2. In this format, /r = —3 corresponds to the point I (0, 
1/3) in Fig. /r = 0 to the Landan point T (1/3,1/9) (LP), and /r = +1 to 
V (0.5, 0.0). In particular, we observe that /i = —1 corresponds to A* ^ 0.57 
located at K (0.2,0.2) in Fig. 

For simulation purposes, the general Hamiltonian in Eqn. ([^ is conve¬ 
niently recast as a biaxial mesogenic lattice model, where particles of D 2 h 
symmetry, represented by unit vectors Vb on lattice sites a and b interact 
through a nearest-neighbour pair potential [15] 


U = -e{G33 - 27 (Gn - ^ 22 ) + A[2(G'n + ^ 22 ) - Gsa]}- (6) 

Here fab= (Va-Ub), Gab=P 2 {fab) with P 2 denoting the second Legendre 
polynomial. The constant e (set to unity in simulations) is a positive quan¬ 
tity setting the reduced temperature = ksT/e, where T is the absolute 
temperature of the system. This is recast along IV of the triangle, using 
Eqn. (17) in reference [10], in terms of the parameter /i as 


H — e[/iGii -h (— 2 G 33 — 2 G 22 + Gii)]. (7) 


3 Details of Simulation 

The Wang-Landau (WL) sampling [12| is a flat histogram technique designed 
to overcome energy barriers encountered, for example, near first order tran¬ 
sitions, by facilitating a uniform random walk along the energy (E) axis 
through an appropriate algorithmic guidance. The sampling, originally de¬ 
veloped for Hamiltonian models involving random walks in discrete configura¬ 
tional space, continues to be applied to various problems in statistical physics 
iaiizi, polymer and protein studies |1H11191|50| and is being developed for 
more robust applications for continuous systems [m E2l [531 El I5511561 EH 
and self assembly [SHI- The proposed algorithm was modified EH] to suit lat¬ 
tice models like the Lebwohl-Lasher interaction [20] , allowing for continuous 
variation of molecular orientations. It was subsequently augmented with the 
so-called frontier sampling technique [Hil to simulate more complex sys¬ 
tems like the biaxial medium. The WL sampling is based on effecting a con¬ 
vergence of an initial distribution over energy E to the density of states (DoS) 
g{E) of the system iteratively. Frontier sampling technique is an algorith¬ 
mic guidance, provided in addition to the WL routine, by which the system 
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is constrained to visit and sample from low entropic regions. The modi¬ 
fied Wang-Landau algorithm using entropic sampling augmented by frontier 
sampling [H] is described below. 

We consider a cubic lattice (size: L x L x L, L = 15, 20) with each lattice 
site representing a biaxial molecule, and hence hosting a (right-handed) triad 
of unit vectors. We initiate the process by assigning random orientations of 
all the axes at every site, and compute the energy of the system at the chosen 
point in the ( 7 , A) plane with the Hamiltonian in Eqn. (|^ (corresponding to 
^ = 1 in Eqn. ([^), under periodic boundary conditions, - the temperature is 
thus measured in reduced units. The energy range of interest of the system 
{Emin, Emax) is divided into N bins (we set N = 40 L^) of equal width, and 
the bin energies are indexed as Ei, corresponding to the values at the centre 
of the bin. We indexed these bins starting from Emin- We initialise g{Ei) 
to an array g^'^\Ei) with equal values {i = 1, where the superscript 

is the iteration run index and the subscript is the energy bin index. The 
estimate of g{Ei) is improved by updating iteratively, until it converges to 
the density of states within a set tolerance limit. 

For liquid crystal systems with continuous degrees of freedom for the ran¬ 
dom walk in conhguration space, we hnd it necessary to perform the simula¬ 
tions on a log-log scale to avoid issues of large numbers and consequent over¬ 
flow problems. Following [61], we work with Q = log(Q:i) = log(log(( 7 (Ej))), 
where ai represents the microcanonical entropy. The acceptance criterion as 
well as reweighting procedures are implemented on this scale. 

During the random walk, the system is permitted to transit from an initial 
conhguration with an instantaneous value Q to a trial conhguration with C,t 
with a probability given by 

p = mm{l,exp[-exp[Ct -h log(l - exp(-(Ci - Cc)))]]}- (8) 

We update the values of Q {i = l,....,iV) of the bins with a Gaussian 
centered at the accepted bin energy value (say Eq), as 

0 ^ Ci+ 7 oexp(—('g) 

Here (70, <5) represent the modihcation parameters. We kept 5 constant 
through the simulation (at 0.002 * N) and chose the initial value of 7 = 0.1. 
Random walk of the system over the energy bins, at this value of 70, is car¬ 
ried out for a large number of lattice sweeps (attempted moves), typically 
10^ sweeps or more depending on the system size. The 70 value is reduced to 
7o —)• 0 . 9570 , and the procedure is repeated until 70 reaches a set small value 
close to zero ( 10“^). The computations involving a progressive reduction of 



7 o, starting from the initial high value, to the set low value constitute an iter¬ 
ation. After two such successive iterations, the differences between histogram 
values at each bin are determined. If the differences are nearly uniform over 
some energy range, it implies that this region is adequately sampled. We 
expect, on entropic grounds, that the flatness of the distribution, in terms 
of fairly uniform increments of histogram values, is achieved more readily 
starting from the maximum energy value. We call the limiting lower energy 
bin, satisfying the flatness criterion, as the frontier, say Ec- Following the 
suggestion of [ 12 ], we update the values of the histogram above Ec by a uni¬ 
form value (say, 0.5). This makes the system, under the above acceptance 
criterion, to perform random walk preferentially in the lower energy region 
hosting less accessible states,- until the histogram values build up to match 
the values in the higher energy regions, above Ec- This process is continued 
with such iterations, and new frontiers are identihed at progressively lower 
energy values, corresponding to an approximate estimation of the DoS over 
larger energy ranges, until the frontier reaches Emin- 

Consequently, a long smoothing run is performed (no frontiers are identi¬ 
hed at this stage) starting with initial values of ( 70 , 5) set to (0.001, 0.002*A^) 
and the value of 70 is progressively decreased during this computation until 
it reaches practically zero value, ~ 10“®. Such iterations continue until a 
specihed hatness criterion is met over the entire energy range. This ensures 
that the hnal C(-^i) converges to its asymptotic value and is representative of 
the density of states of the system, within the tolerances prescribed by the 
hatness criterion. 

We now construct a large entropic ensemble of microstates (say, M ~ 
4 X 10^) by ehecting a random walk of the system over the energy bins 
(f = 1,...,A^) with an acceptance probability based on (~^{Ei) (analogous 

to Eqn. (|^). We label the microstates as C* [i = = 1,.M] 

with M ^ N. We note that an bin for example hosts a large number of 
microstates (C*) with distinct energies E{CD, however represented by the 
same density of states g^Ef). 

The relevant thermodynamic quantities are calculated at each tempera¬ 
ture by constructing appropriate canonical ensemble of states using a reweight¬ 
ing technique [221 • We refer to these ensembles as RW-ensembles, to diher- 
entiate from those constructed through the Metropolis guided random walk 
(B-ensembles). The equilibrium averages of a physical variable ’O’ at a tem¬ 
perature T {1^=-^) are computed through this procedure as 


{O) 


Ec.Q(Q)i7(i^.)exp [-mCl)] 
c/(Ei)exp[-/?E(C'*)] 


( 10 ) 
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The representative free energy F, as a function of the energy of the system, 
as well as of the two dominant order parameters (uniaxial and biaxial orders) 
is computed from the DoS and the microcanonical energy, - both available 
as a function of bin number in the entropic ensemble. 

The WL simulations were carried out at different values of ( 7 , A) in 
Eqn. ([^ so as to trace the trajectory OIV of the essential triangle in Fig. 
at about 60 chosen points. For purposes of comparison, conventional MC 
sampling (based on Metropolis algorithm) was used to construct canonical 
(Boltzmann) ensembles. Considering an attempted N= moves as one lat¬ 
tice sweep (MC step), the system is equilibrated, and a production run is 
carried out, each for 6 x 10® MC steps. In our analysis, we hud it necessary 
to distinguish between the averages from B-ensembles and RW-ensembles. 

The physical parameters of interest in this system, calculated at each A*, 
are the average energy < E >, specihc heat < >, energy cumulant V4 

(= 1— < > /(3 < >^)) which is a measure of the kurtosis [63], the 

four order parameters of the phase calculated according to [niiMi and their 
susceptibilities. These are the uniaxial order < Fqq > (along the primary 
director), the phase biaxiality < molecular contribution to 

the biaxiality of the medium < R 22 >, and the contribution to uniaxial order 
from the molecular minor axes < Rq 2 >■ 

The averages are computed at a temperature resolution of 0.002 units 
in the temperature range [0.05, 2.05]. The temperature of the simula¬ 
tion is scaled to conform to the values used in the mean held treatment: 
1/P* = 37'/(9[2F(l + 3A)]) [71 Ej. Statistical errors in different observables 
are estimated over ensembles comprising a minimum of 5 x 10® microstates, 
and these are compared with several such equilibrium ensembles at the same 
( 7 , A) value, but initiating the random walk from different arbitrary points 
in the conhguration space. We hnd the relative errors in energies are 1 in 
10®, while those in the estimation of the order parameters are 1 in 10^. We 
also note that these error estimates from RW- ensembles are smaller relative 
to B-ensembles of comparable size by at least an order of magnitude, ow¬ 
ing to the efficacy of the importance sampling involved in the reweighting 
procedure. 


Results 


The temperature variations of the specihc heat and the two dominant scalar 
order parameter (Fqq and RI 2 ) values obtained from RW-ensembles at various 


values of A* along the arc OIT (A*-axis) are shown in Figs. 2(a) - 2(d) 


It is noted from Fig. 2(a) that for all values of A* in the range 0.1 
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(c) (d) 

Figure 2: (color online) Temperature variation of the specific heat (in arbi¬ 
trary units) for different ranges of A*: (a) 0.1 - 0.25; (b) 0.26 - 0.33; (c) 0.33 
- 0.53; (d) 0.53 - 0.733. Corresponding variations of the two primary order 
parameters {Rqq and 7 ^ 22 ) are shown in the insets with same colour scheme. 
Splitting of the Cy peaks in (d) for AJ^> 0.53 and qualitative changes in the 
temperature variation of order parameters are clearly observed. 






























































0.25, two transition peaks are observed in the specific heat. As the biaxial 
system is cooled from the high temperatnre isotropic phase, an initial I — Nu 
transition occnrs at a high temperatnre Ti followed by a second transition 
Nu—Nb at lower temperatnre T 2 . The I—Nu transition temperatnre remains 
fairly constant with the variation in A*, whereas Nu — Nb transition shifts 
towards higher temperatnres as A* increases from 0.1 to 0.25. This behavionr 
is also reflected in the order parameter profiles shown in the inset. The two 
transitions eventually coalesce at A* = 0.26 resulting in a triple point and a 
direct isotropic-biaxial (/ — Nb) transition occurs from A* = 0.26 to 0.53, as 
depicted by the specific heat profiles and order parameters (inset) of Figs. 
|2(b) - [2(^ These results from RW-ensembles agree qualitatively with those 
obtained from the B-ensembles in the range of A* = 0.1 — 0.53. 

A comparative study of the WL and MC simulation results for certain 
representative values of A* are shown in Figs. 3(a) - 3(d) It is observed that 


qualitative agreement with the mean-field predictions exists upto A* < 0.53 
and deviations of the RW-ensembles from MF and B-ensembles start from 
A* = 0.54(5/29,19/87) (i.e C 3 in the essential triangle of FigQ. 

the results from RW-ensembles agree with MF 


Referring to Fig. 2(a) 


predictions except for the actual values of the location of the tricritical and 
triple points Ci and C 2 - In this respect, one has to make allowances for 
unavoidable finite size effects on the simulation data on the one hand, and 
the inherent approximate nature of the mean field theoretical analysis in this 
respect, on the other. 

At L=20, the simulation results show that the tricritical point lies in the 
neighbourhood of A* = 0.18; the nature of the Nu — Nb transition appears 
to change to a (weak) first order for values of A* > 0.18, as evidenced from 
the energy cumulant data shown in Fig. The triple point is located at 
A* ~ 0.26 (corresponding MF value is ~ 0.22) as the transition sequence 
/ — Nu — Nb changes to / — Nb at this value of A* (see Fig. 2(b)). Transition 
temperatures derived from these simulations are summarised in Table 

It is of interest to observe that the I—Nb transition progressively becomes 
very strong first order as A* value increases to 1/3 and is most pronounced 
at the point corresponding to coordinates (0,1/3) (Fig. [^. Fig. [^depicts the 
specific heat with (inset) energy cumulant and order parameters with (inset) 
their susceptibilities for this value of A*. This feature of the transition is also 
demonstrated by the variation of the representative free energy obtained from 
the DoS and bin energies, as a function of the two order parameters. As an 
example, we depict its variations across the transition (A* = 0.330, L=15), in 
Fig.g Observation of such a strong free energy barrier (see the coexistence 
region at T = is supportive of the MF prediction at I. Very similar 

plots result as a function of Rqo 
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(c) 


(d) 


Figure 3: (color online) Comparison of the results obtained from B-ensembles 
(hollow red circles) and the RW-ensembles (hollow black squares): Temper¬ 
ature variation of the specihc heat (in arbitrary units) and corresponding 
variations of the two primary order parameters {R^q and R 22 ) are shown for 
different values of A* in regions Of and IV: (a) 0.2; (b) 0.33; (c) 0.51; (d) 
0.54. The overlap of the corresponding curves (a to c) clearly indicates the 
agreement between the two ensembles upto A* = 0.53 as mentioned in the 
text. Fig. |3(d) shows the qualitative disagreement hrst noticed at A* = 0.54. 
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Figure 4: (color online) Variation of energy cumulant with temperature, at 
different A* values: (a) 0.18; (b) 0.2; (c) 0.22 (d) 0.26 


Table 1: Transition temperatures in the range of A*= (0.18, 0.26): and 

7^2 are transition temperatures (in reduced units) obtained from simulation, 
while T[ and are the corresponding equivalent mean held temperatures. 


A* 





0.18 

1.1753 

0.9919 

0.1272 

0.1074 

0.2 

1.1937 

1.0770 

0.1243 

0.1122 

0.22 

1.2110 

1.1490 

0.1216 

0.1153 

0.26 

1.2516 


0.1172 
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Figure 5: (color online) (a) Specific heat profile with (inset) energy cumulant 
(b) Order parameters with (inset) susceptibility prohles for A*=0.33. Point I 
of the essential triangle is the intersection point of the three uniaxial torque 
axes [9] and hosts the strongest hrst order I — Nb transition, as shown by 
the very sharp features of these physical properties. 
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Figure 6 : (color online) Representative free energy (in arbitrary units) as a 
function of biaxial order parameter R 22 (^) ^ (b) T = Tj^^i and 

T < T^gi at A*=0.33 for a system of size L=15. 


Referring to Fig. 2(d), we observe that the peak splits starting from 


A* ~ 0.54, signalling the onset of two transitions. The temperature gap be¬ 
tween transition peaks increases with A* above this value, attaining a max¬ 
imum at A* = 0.733 (T on the triangle). These observations are illustrated 
in Fig. plotting all the relevant variables as a function of temperature at 
chosen values of A* (0.54, 0.58, 0.62, 0.66, 0.69, 0.72) along the segment C 3 T. 
These graphs depict the temperature variation of specific heat Cy with en¬ 
ergy cumulant V 4 (inset) and order parameter (i?oo> ^ 22 ) profiles along with 
respective susceptibilities y (inset). 

The nature of the two phases below the clearing point is inferred from 
the order parameter prohles and their susceptibility peaks. Referring to 
the two transition temperatures in decreasing order as Ti and T 2 , the data 
indicate that the onset of a biaxial phase takes place at Ti itself, and the 
growth of biaxial order in the intermediate phase is marginal as compared to 
the uniaxial order. Further, both the uniaxial and biaxial order parameters 
display a sudden upward jump at T 2 , and subsequently increase rapidly (more 
pointedly the biaxial order R 22 ) ^^e temperature is lowered further. This 

behaviour is prominent in the neighbourhood of A* = 0.66. The susceptibility 
of i?QQ exhibits two peaks corresponding to the two transitions, whereas that 
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(a) 


(b) 



(c) (d) 




(e) (f) 
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Figure 7: (color online) Specific heat profile with (inset) energy cumulant 
V 4 and order parameters with (inset) susceptibility profiles at different A* 
values (a) 0.54; (b) 0.58; (c) 0.62; (d) 0.66; (e) 0.69; and (f) 0.72 










































(a) 


(b) 


Figure 8 : (color online) Free energy shown as a function of (a) and (b) 
i? 22 ) the transition temperature Ti for A* values in the region C 3 T of 

Fig.[g 


of RI 2 shows only a single peak at T 2 , for all values of A*. The energy 
cumulant V 4 shown in the inset of each of the figures indicates the first order 
nature of the I — Nbi transition (at Ti). The additional second dip at the 
lower temperature transition (at T 2 ) appears to point towards the progression 
of the hrst order nature of the Nbi — Nb transition. It is observed that the 
dip in the cumulant at T 2 is maximum at A*= 0.66. 

We also examined the representative free energy plotted as a function 
of the order parameters at the transition temperatures Ti and T 2 . These 
variations observed near Ti, for different A* values (covering the region CsT) 
are shown in Figs. Focussing on Figs. 8 (a) and 8 (b), one immediately 
observes that both the free energy profiles (with respect to Rqq and R 22 ) at Ti 
show a distinctive indication of a developing minimum at a lower temperature 
evidenced by the systematic deviations (from a smooth continuation) of the 
profiles at the respective higher values of the two order parameters. And, 
as A* increases in the region, the location of these sharp deviations 
progressively shift to a higher value of the corresponding order parameter. 

We tracked the variation of these profiles closely from Ti to T 2 (shown 
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(a) (b) 

Figure 9: (color online) Free energy shown as a function of (a) i?oo and (b) 
oil cooling from Ti to T 2 for A* = 0.65 
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(a) 


(b) 


Figure 10: (color online) Free energy shown as a function of (a) Rq^ and 
(b) R 22 , at the transition temperature T 2 for A* values in the region (FsT of 


Fig. [3 
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in Fig. I^for a single value of A*=0.65), and find that the free energy minima 
gradually shift towards high order regions, and the second transition at T 2 
corresponds to a gradual shift of the free energy minima towards the curious 
regions, depicted in Fig. The variation of the free energy at T 2 is shown in 
Figs. 10(a) I and 10(b) for different values of A* (in region C 3 T). These depict 
a free energy minimum attained at T 2 for all values of A*. 

We argue that the progression of the free energy prohles with temperature, 
as a function of and RI 2 , and matching of the values of respective order 
parameters at T 2 with the location of sharp deviations observed in Fig. are 
further evidences for the existence of two transitions in this region. It may 
be pointed this could be made possible only by adopting a MC sampling 
procedure which facilitates the computation of free energy prohles of the 
system, via the density of states. 

We are thus led to the conclusion that in this region of A* values, the 
medium undergoes two transitions, and both the low temperature phases 
have biaxial symmetry. From the data on the limited temperature region 
available for the intermediate phase, and in comparison with the low tem¬ 
perature phase, it appears that the biaxial order in the intermediate phase is 
somewhat inhibited, presumably by free energy barriers. It is only after the 
second transition at T 2 (between the two biaxial phases, and hence necessar¬ 
ily a hrst order transition) that the biaxial order shows normal increase with 
decrease in temperature, as is to be expected. Thus we propose the phase 
sequence in this region to be Nb — Nbi — I- 

We now construct the phase diagram as a function of the arc length A* 
based on the specihc heat data (from RW-ensembles), shown in Fig. 11 at 56 
values of A* distributed over the arc 01V (see [39] for details). The transition 
temperatures at a few representative values of A* beyond the Landau point 
T (segment TV) are obtained from the B-ensemble data. The temperature 
of the simulation is scaled to conform to the values 1/(3* used in the mean 
held treatment as discussed in section III. 

A comparison of the phase diagram proposed from the current MC sim¬ 
ulations [39| with the one predicted based on mean-held theory [9] brings 
out clear qualitative diherences in the region CaTV of the essential triangle. 
We observe that the predicted direct transition from the isotropic to biaxial 
phase is replaced by two transitions in which an intermediate biaxial phase 
occurs between these two phases. These results start deviating starting from 
A* > 0.54, very close to the point C 3 in Fig. The fact that B-ensembles 
constructed from simple conhgurational random walks based on Metropolis 
algorithm fail to detect the second transition in this A* region merits some 
discussion. 
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Figure 11: (color online) Phase diagram as a function of A*, derived from RW- 
ensembles. The transition temperature 1/(3* is scaled to conform to mean - 
held values as indicated in the text. Points along OIV in Fig. [^are mapped 
onto the A*-axis for reference. An additional biaxial-biaxial transition is 
observed in the region KTV in place of a single transition (to the biaxial 
phase) predicted by the mean-held theory [3^ . 
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Figure 12: (color online) Comparison of data, as a function of temperature, 
from the B- and RW- ensembles, at the Landau point T (1/3,1/9): (a) 
Specihc heat and (b) order parameters Rqq and -R 22 - The insets focus on 
(a) the energy cumulant 14 and (b) order parameter susceptibilities (y's), 
both derived from RW-ensembles (A* = 0.733). 


Discussion 


In order to look for the origin of the additional low temperature specihc heat 
peak observed in the region C 3 T, which was not detected by Boltzmann sam¬ 
pling, we made a comparison of the simulation results from RW-ensembles 
with those obtained from B-ensembles at A*=0.733 (1/3,1/9) (Landau point 
T) shown in Fig. 12 The location of T is unique as it is the intersection 


point of the dispersion parabola with the segment IV. MF theory predicts a 
direct transition from the isotropic to biaxial phase at this point, and it is 
the only such point on the parabola. From the perspective of the interac¬ 
tion Hamiltonian, the coordinates ( 7 , A) of T represent a unique symmetry: 
The Hamiltonian has no interaction between the uniaxial components of the 
molecular tensor (/x = 0 at T, see Eqn. ([^), and is purely biaxial in nature 
(involving m and axes). A consequence of the present Endings in this 
context is curious. They do confirm the presence of a direct transition from 
the isotropic to biaxial symmetry, but these also indicate that there is yet 
another biaxial-to-biaxial transition at a lower temperature. Further, the 
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onset of the first biaxial phase at Ti is not leading to a natnral progression 
of the biaxial order with decrease in temperatnre, and it is only after the 
transition at T 2 that the macroscopically signihcant and hence observable, 
R 22 valne seems to be realizable. 

We thns focns on the Landan point, and present the simnlation resnlts ob¬ 
tained from the two types of MC sampling methods: data from B-ensembles 
and from RW-ensembles. Fig. 12 shows the specihc heat (energy cnmulant as 
inset) and order parameters (snsceptibilities as inset), compnted as a fnnc- 
tion of temperatnre at the point T, obtained from these ensembles. It may be 
noted that the derived physical variables from B-ensembles do not betray the 
onset of the second transition at T 2 , thus lending support to MF predictions, 
as has been noted in the earlier report on this work |39j . 

We now examine the contour maps of the distribution of microstates in 
the entropic ensemble (set of microstates which are approximately uniformly 
distributed with respect to energy) collected at the Landau point. Fig. [T^a) 
depicts such a contour map in the space of uniaxial order and energy (per 
site), along with the thermal averages computed from RW-ensembles and B- 
ensembles superposed for ready comparison. Similarly, Fig. 13 b) shows cor¬ 
responding contour map plotted between biaxaial order and energy (per site), 
along with thermal averages of the two canonical ensembles again superposed. 
The traversal path of the B-ensemble averages is seen to be encompassing 
regions corresponding to contour peak positions, whereas the RW-ensemble 
average is observed to follow a different trajectory, consequent to encompass¬ 
ing a wider collection of microstates visiting sparse regions, corresponding to 
large deviations of the order parameter. This is seen as a manifestation of the 
process of collection of microstates of the entropic ensemble by the algorithm 
employed, representative in their distribution (with respect to energy) of the 
underlying density of states. As has been pointed out and argued earlier 
(Fig. 6 in [3^]), the algorithmic guidance of the WL-procedure is seeking out 
all accessible microstates (an approximate microcanonical ensemble) in each 
bin of energies, in the process visiting relatively rare states which correspond 
to larger excursions in the order parameter, and hence correspondingly larger 
fluctuations of the component energies of the Hamiltonian in Eqn. ([^, while 
conforming to the same energy bin. The requirement of an accurate determi¬ 
nation of DoS through the entropic sampling procedure apparently demands 
inclusion of these rare microstates, and the process of reweighting used to 
construct the equilibrium ensembles through this elaborate procedure, in¬ 
cludes them in the thermal averages as a consequence. 

Thus the differences observed in the averages from the two procedures 
are to be appreciated from the standpoint of simulations. As has been dis¬ 
cussed |39] , these rare microstates indeed correspond to situations where the 


24 





-3 -2 -1 


E 


Figure 13: (color online) Contour plots of the distribution of microstates 
collected in the entropic ensemble at A* ~ 0.733: (a) Microstate energy 
versus its uniaxial order and (b) Microstate energy versus its biaxial order. 
The superimposed red (dash dotted line) and black (dashed) lines are thermal 
averages from B- and RW-ensembles, respectively. 
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ordering of either of the molecular axes (involved in the D^h symmetry of the 
pair-wise interaction, i.e m or ej_ ) form a spontaneous and equally probable 
calamitic axis during the evolution of the system, by virtue of having the 
largest instantaneous eigenvalue of the diagonalised ordering tensors of the 
three molecular axes. It is apparent that conventional sampling methods, 
not under algorithmic compulsion to estimate the DoS of the system, are not 
geared to sample such rare states, and hence come up with different aver¬ 
ages. In the process, it appears that second low temperature transition is 
not evident in the earlier work. 

In addition, we note that the deviation of the simulation results from 
the mean held expectations occur along the diagonal IV of the essential tri¬ 
angle where the pairwise interaction Hamiltonian has symmetry and 
is expressed in the reduced form as in Eqn. (|^ in terms of a single pa¬ 
rameter fi. It is observed that the deviations start from A* > 0.54 (point 
Cs (5/29,19/87) and continue till the Landau point where /i = 0. It may be 
noted that the point K (0.2, 0.2) corresponding to /i = —1 is very close to 
the point C 3 (0.172,0.218). It can be inferred from Eqn. ([^ that, starting 
from the neighbourhood of K, the uniaxial attractive coupling of the e-axis 
becomes lower in strength than that of the (biaxial) attractive coupling of 
the other two axes, and continues to decrease as A* increases on the diagonal. 
As a result, the ordering of the biaxially coupled e± and m axes is favoured 
as temperature is decreased, leading to the hrst onset of a biaxial symmetry 
from the isotropic phase. This is followed by an ordering of the e-axes at a 
lower temperature, leading to the stabilization of both the orders, in partic¬ 
ular biaxial order. Thus it appears that the growth of biaxial order in the 
intermediate biaxial phase is inhibited by the lack of long range order of the 
e-axes (in this region of A*-axis). As one approaches the Landau point, /i 
in Eqn. ([^ tends to zero (from the negative side), thereby suppressing the 
second transition temperature as well as weakening the efficacy of this term 
to drive a transition. This lack of concomitant ordering of all the molecular 
axes leads to inhomogeneity in the medium, and we tend to attribute all the 
interesting aspects of the simulation to this feature of the Hamiltonian. 

Finally, we wish to comment on the curious role played by the WL- 
algorithm in the analysis of the phase diagram. It has been already estab¬ 
lished that this algorithm assists the system in overcoming energy barriers 
of the system, as the simulation pushes the system to make a random walk 
in the configuration space which is uniform with respect to energy. A suc¬ 
cessful convergence of the probability density yields a limiting distribution of 
microstates with respect to the total energy of the system - the representative 
density of states. The role of this algorithm in the present study seems to be 
qualitatively different and yet illustrative of its varied applicability. The WL 
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algorithm, even while operating within a single energy bin (an approximate 
microcanonical ensemble), appears to seek ont rare states, corresponding to 
the otherwise inaccessible flnctnations of the component energies, making up 
the total energy (see Fig. 6 in [39]). Inclusion of these microstates ( as 
representative states for purposes of computing averages) is naturally em¬ 
bedded in the WL-method, while estimating the DoS accurately. We argue 
that the Metropolis sampling fails to access these states due to apparent en¬ 
ergy barriers within the system inhibiting sampling of microstates with such 
large fluctuations in their energy components. We conclude that the new re¬ 
sults reported here are the outcome of this facet of efficiency of the entropic 
sampling. 

6 Conclusions 

In conclusion, we present compelling evidences from Monte Carlo simulations 
based on entropic sampling, to propose an additional biaxial phase along a 
region of the arc of the essential triangle, augmenting our earlier report [39] . 
The arguments advanced in this respect, particularly of the inevitable pres¬ 
ence of inhomogeneities in the absence of a long-range order of the third 
stabilising axis e, seem to lend support to the findings (based on Boltzmann 
MC sampling) reported in the partly repulsive region of the A*-axis (seg¬ 
ment TV: the /i-model [IQ]). At a more general level we conclude that the 
cross-coupling between the uniaxial and biaxial tensor components of the 
neighbouring molecules (y-term in Eqn. ([^) seems to be playing an impor¬ 
tant role in determining the phase sequences. Further, we suggest that its 
signihcant presence, even along trajectories inside the triangle (which could 
be relevant for practical purposes) should have such inhibitive influence on 
the condensation of a biaxial phase with measurable biaxial order. Our re¬ 
cent simulational work on two such trajectories interior to the triangle are 
supportive of this conjuncture (to be published). 
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